Blight Predictions for Detroit

In this project, I will attempt to predict whether, for any given parcel in Detroit with a building on it, whether a building on that parcel will be be targetted for demoltion. Potential predictors include citations relatated to the building, crime, and complaints concerning the building related to blight. The probject is currently at the data cleaning stage, and I may add some other data to the project. So far I am using about 3GB of data downloaded from https://data.detroitmi.gov/.

library(tidyverse)
library(sf)
library(ggmap)
library(lwgeom)
#recorded violations associated with blight (e.g. unkempt properties)
blight_violations <- read_csv("./data/Blight_Violations_3_19_2018.csv", 
                              guess_max = 10^6)
#read the downloaded file for all the building permits and then filter out the permits for dismantling
dismantle_permits <- read_csv("./data/Building_Permits_3_19_2018.csv", 
                             guess_max = 10^6) %>%
  filter(`Building Permit Type` == "Dismantle")
  
#the files that contain the crime data
crime_to_12062016 <- 
  read_csv("./data/DPD__All_Crime_Incidents__January_1__2009_-_December_6__2016.csv",
           guess_max = 10^6)
crime_12062016_to_03192018 <- 
  read_csv("./data/DPD__All_Crime_Incidents__December_6__2016_-_3_19_2018.csv", 
           guess_max = 10^6)
#the 311 system
improve_detroit_issues <- read_csv("./data/Improve_Detroit_Issues_3_19_2018.csv", 
                                   guess_max = 10^6)
#another file with demolition information downloaded 4/4/2017
completed_demolitions <- read_csv("./data/Detroit_Demolitions.csv",
                                  guess_max = 10^6)
#the shapefile representing Detroit parcels, read into 
parcel_sf <- st_read("./data/Parcel Map")
#convert the parcel number column to a character vector, to match `Parcel Number` in the dismantle permits data
parcel_sf <- parcel_sf %>% mutate(parcelnum = as.character(parcelnum))

For all of the downloaded datasets other than the parcels dataset, we extract the usable latutude and longitude values and then use this information to form simple features (sf) objects. Rows with obviously incorrect values, or values that would represent positions well outside Detroit, are filtered out, together with rows for which the latitude or longitude data is missing.

#function for converting the position (character) column into a column of points in the simple features (sf) framework.
add_sf_point <- function(df, column) {
  
  #extract the latitude and longitude from the string column that contains both. With the parentheses
  #located from the end of the strings, it is possible to use the the same function for all five of 
  #the datasets for which we need to extract this information.
  latitude <- str_sub(df[[column]], 
                      stringi::stri_locate_last_fixed(df[[column]], "(")[,2] + 1,
                      stringi::stri_locate_last_fixed(df[[column]], ",")[,1] - 1)
  longitude <- str_sub(df[[column]], 
                       stringi::stri_locate_last_fixed(df[[column]], ", ")[,2] + 1, 
                       stringi::stri_locate_last_fixed(df[[column]], ")")[,1] - 1)
  
  #add the latititude and and longitude to a copy of the dataframe, filter out the NAs from
  #these results, and then convert convert to sf, with point positions indicated in the
  #geometry column
  mutated <- df %>% mutate(extracted_lat = as.double(latitude),
                      extracted_lon = as.double(longitude))
  
  #remove rows with NAs for latitude or longitude, or with values well outside of Detroit
  filtered <- mutated %>%
    filter(!is.na(extracted_lat) & !is.na(extracted_lon)) %>%
    filter(41 < extracted_lat & extracted_lat < 44 & -85 < extracted_lon & extracted_lon < -81)
    
  #create a dataframe from the items that have been filtered out
  result_coord_na <- setdiff(mutated, filtered)
  
  #create sf objects from the rows with usable latitude and longitude information
  result_sf <- st_as_sf(filtered, coords = c("extracted_lon", "extracted_lat"), crs = 4326)
    
  return(list(result_sf, result_coord_na))
}
#apply the function to the five datasets for which the data was not loaded as a simple features dataframe, thus producing a list of two dataframes for each of the datasets, the first element of the list a simple features data frame and the second element a dataframe with the instances for which it was not possible to convert to simple features
blight_violations_split <- add_sf_point(blight_violations, "Violation Location")
dismantle_permits_split <- add_sf_point(dismantle_permits, "Permit Location")
crime_to_12062016_split <- add_sf_point(crime_to_12062016, "LOCATION")
crime_12062016_to_03192018_split <- add_sf_point(crime_12062016_to_03192018, "Location")
improve_detroit_issues_split <- add_sf_point(improve_detroit_issues, "Location")
completed_demolitions_split <- add_sf_point(completed_demolitions, "Location")

We now consider the data for which we do not yet have position data, and complete the information as well as we reasonably can, using the Google api and a function, geocode_pause, that handles some of api’s quirks.

#the portino of the downloaded blight citations data, for which we do not have
blight_vio_na <- blight_violations_split[[2]]
#remove the rows for which geocoding is not likely to prodoce reliable results
useful <- blight_vio_na %>% filter(!is.na(`Violation Street Name`), 
                                   `Violation Street Number` > 0, 
                                   !is.na(`Violation Zip Code`))
#create a column of addresses to be used in geocoding
useful <- useful %>% 
  mutate(complete_address = paste(`Violation Street Number`, " ", `Violation Street Name`, ", ",
                                  "Detroit, Michigan", " ", `Violation Zip Code`, sep = ""))
#function makes a maximum 6 attempts to geocode the given address using the Google API, with a pause of 1 second between attempts. We will use the function for the other datasets as well.
geocode_pause <- function(address) {
  for (index in 1:6) {
    Sys.sleep(1)
    location <- ggmap::geocode(address)
    if (!is.na(location$lon)) {
      return(location)
    }
  }
}

#apply geocode_pause to each of the elements of the complete_addresse column and place the result in a new column, in which each entry is a data frame 
useful <- useful %>% mutate(location = map(complete_address, geocode_pause))

#save to disc, to avoid avoid the need to geocode these addresses again when we rerun the analysis
write_rds(useful, "./data/blight_violations_geocodes.rds")

The geocoding has returned a data frame for each of the addresses. We thus need to unpack the elements of the location column, each of which is a data frame.

#read blight_violations_geocodes as a tibble
blight_violations_geocodes <- read_rds("./data/blight_violations_geocodes.rds")
#function for removing the instances for which geocoding failed (for which the value in the location column is NULLL). We will use this function for all of the geocoded data frames.
remove_null_locations <- function(df) {
  #identify the rows for which the value in the location column is NULL
  null_rows <- list()
  for (index in 1:nrow(df)) {
    if (is.null(df$location[[index]])) {
      null_rows <- c(null_rows, index)
    }
  }
  #remove the rows for which the value of the location column is NULL
  df <- df[-as.integer(null_rows),]
}
blight_violations_geocodes <- remove_null_locations(blight_violations_geocodes)
#With blight_violations_geocodes a tibble, we can apply tidyr::unnest(), which will place the latitude and longitude in columns labelled "lat" and "lon".
blight_violations_geocodes <- blight_violations_geocodes %>% unnest(location)
#fill in the `Violation Latitute` and `Violation Longitude` data frames, which alread exist in the blight_violations data frame
blight_violations_geocodes <- blight_violations_geocodes %>%
  mutate(`Violation Latitude` = lat,
         `Violation Longitude` = lon)
#cut out some columns that have been added
blight_violations_geocodes <- blight_violations_geocodes %>% 
  select(-extracted_lat, -extracted_lon, -complete_address)
#put the position information into a simple features format (which will remove the "lat" and "lon" columns)
blight_violations_geocodes_sf <- st_as_sf(blight_violations_geocodes, 
                                          coords = c("lon", "lat"),
                                          crs = 4326)
#combine the results with the previously generated sf data
blight_violations_sf <- rbind(blight_violations_split[[1]], blight_violations_geocodes_sf)
rm(blight_vio_na, blight_violations, blight_violations_geocodes, 
   blight_violations_geocodes_sf, blight_violations_split, useful)

#the dismantle permits for which position data (latitude and longitude) is missing
dismantle_permits_split_na <- dismantle_permits_split[[2]]
#remove the last two columns, which were not contained in the original dismantle_permits datastet
dismantle_permits_split_na <- dismantle_permits_split_na %>% 
  select(-extracted_lat, -extracted_lon)

#geocode the items in dismantle_permits_split_na, using the address column and the function geocode_pause, which makes a maximum of six attempts for each item. The result is list of dataframes in the location column.
dismantle_permits_split_geocode <- dismantle_permits_split_na %>%
  mutate(location = map(str_c(`Site Address`, ", Detroit, Michigan"), geocode_pause))

#write the results of the geocoding to disk, to avoid having to repeat the geocoding when rerunning the analysis.
write_rds(dismantle_permits_split_geocode, "./data/dismantle_permits_geocodes.rds")

rm(dismantle_permits_split_na)

#load the geocoded data frame into R
dismantle_permits_split_geocode <- read_rds("./data/dismantle_permits_geocodes.rds")
#use the remove_null_locations() to remove the rows for which geocoding failed and then parse the information  in the dataframes in the location column into two new columns, lat and lan
dismantle_permits_split_geocode <- 
  remove_null_locations(dismantle_permits_split_geocode) %>%
  unnest(location)
#convert to a simple features (sf) data frame, using the latititudes and longitudes
dismantle_permits_geocode_sf <- st_as_sf(dismantle_permits_split_geocode, 
                                          coords = c("lon", "lat"),
                                          crs = 4326)
#append this simple features dataframe to the dataframe for which we already had usable positions
dismantle_permits_sf <- rbind(dismantle_permits_split[[1]], dismantle_permits_geocode_sf)
rm(dismantle_permits_split_geocode, dismantle_permits_geocode_sf, dismantle_permits, dismantle_permits_split, dismantle_permits_split_na)

We now fill-in the missing position information for the dataset for crimes up to 12-06-2016

#return to the older crime data
crime_to_12062016_leftovers <- crime_to_12062016_split[[2]]
#cut out the addresses that begin with "00"
crime_to_12062016_leftovers <- crime_to_12062016_leftovers %>% 
  filter(str_sub(LOCATION, 1, 2) != "00")
#filter out some obviously useless addresses, with few characters before the first "("
crime_to_12062016_leftovers <- crime_to_12062016_leftovers %>% 
  filter(!str_locate(LOCATION, "\\(")[,1] %in% 1:13)
#remove the two columns that were added earlier
crime_to_12062016_leftovers <- crime_to_12062016_leftovers %>% 
  select(-extracted_lat, -extracted_lon)
#create a column for use in geocoding
crime_to_12062016_leftovers <- crime_to_12062016_leftovers %>% 
  mutate(extracted_address = str_c(str_sub(LOCATION, 1, 
                                           str_locate(LOCATION, "\\(")[,1] - 2),
                                                      ", Detroit, Michigan"))

#geocode the elements of extracted_address, using the function geocode_pause
crime_to_12062016_leftovers_geocode <- crime_to_12062016_leftovers %>% 
  mutate(location = map(extracted_address, geocode_pause))

#save the results, to avoid having to geocode again when rerunning the analysis
write_rds(crime_to_12062016_leftovers_geocode, "./data/crime_to_12062016_leftovers_geocode.rds")

crime_to_12062016_leftovers_geocode <- read_rds("./data/crime_to_12062016_leftovers_geocode.rds")
#cut out the column we used for geocoding
crime_to_12062016_leftovers_geocode <- 
  crime_to_12062016_leftovers_geocode %>% select(-extracted_address)
  
#cut out of the geocode failures and put the location information into the columns lat and lon
crime_to_12062016_leftovers_geocode <- 
  remove_null_locations(crime_to_12062016_leftovers_geocode) %>%
  unnest(location)
#create a simple features (sf) object, using the latititudes and longitudes
crime_to_12062016_leftovers_sf <- st_as_sf(crime_to_12062016_leftovers_geocode, 
                                          coords = c("lon", "lat"),
                                          crs = 4326)
#append this simple features dataframe to the dataframe for which we already had locations
crime_to_12062016_sf <- rbind(crime_to_12062016_split[[1]], crime_to_12062016_leftovers_sf)
rm(crime_to_12062016, crime_to_12062016_leftovers_geocode, crime_to_12062016_leftovers_sf, crime_to_12062016_leftovers, crime_to_12062016_split)

#consider the examples in the recent crime data for which the conversion to sf didn't work, remove the two columns that we have added, and create and address column for geocoding
crime_12062016_to_03192018_leftovers <- crime_12062016_to_03192018_split[[2]] %>%
  select(-extracted_lat, -extracted_lon) %>%
  mutate(extracted_address = str_c(`Incident Address`, ", Detroit, Michigan"))

crime_12062016_to_03192018_geocode <- crime_12062016_to_03192018_leftovers %>% 
  mutate(location = map(extracted_address, geocode_pause))

write_rds(crime_12062016_to_03192018_geocode, "./data/crime_12062016_to_03192018_geocode.rds")

crime_12062016_to_03192018_geocode <- read_rds("./data/crime_12062016_to_03192018_geocode.rds") %>%
  select(-extracted_address)
#remove the rows for which the value of location is NULL and then unnest the remaining locations
crime_12062016_to_03192018_geocode <- 
  remove_null_locations(crime_12062016_to_03192018_geocode) %>%
  unnest(location)
#convert the dataframe to a simple features set
crime_12062016_to_03192018_sf <- st_as_sf(crime_12062016_to_03192018_geocode, 
                                          coords = c("lon", "lat"),
                                          crs = 4326)  
#combine the geocoded data with the sf dataframe created earlier
crime_12062016_to_03192018 <- rbind(crime_12062016_to_03192018_split[[1]], crime_12062016_to_03192018_sf)
rm(crime_12062016_to_03192018_split, crime_12062016_to_03192018_geocode, crime_12062016_to_03192018_leftovers, crime_12062016_to_03192018_sf)

#geocode the one item in the Improve Detroit Issues data for which the given coordinates were obviously incorrect, and then convert to an sf object. If geocoding fails, run this bit again
improve_detroit_issues_leftover_sf <- improve_detroit_issues_split[[2]] %>%
  select(-extracted_lat, -extracted_lon) %>%
  mutate(location = map(Address, geocode_pause)) %>%
  unnest(location) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)  

Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=14530%20%20Vaughan%20Detroit,%20Michigan&sensor=false

#splice with the previously generated sf dataframe
improve_detroit_issues <- rbind(improve_detroit_issues_split[[1]], improve_detroit_issues_leftover_sf)
rm(improve_detroit_issues_split, improve_detroit_issues_leftover_sf)

#create the other set of demolition information
completed_demolitions_sf <- completed_demolitions_split[[1]]
#note that location information in this dataset is complete
completed_demolitions_split[[2]]

We begin the assignment of labels to the buidings: blighted or not blighted. Buildings will be represented by parcels that have or have had buildings on them, whether by being so represented as in the parcels_sf data frame as including structures or in the dismantle permits dataframe as having had a dismante permit associated with it, thus suggesting that there was a building on the parcel.

We will use parcel numbers to refer to the parcels. However, as the following bit of code shows, the parcels dataset contains a few rows in which the parcel humbers are the same (duplicate_parcel_numbers_in_parcel_data contains 78 rows).

parcel_sf <- parcel_sf %>%
  mutate(row_num = row_number())
#As per above, following returns a 78-row data frame
duplicate_parcel_numbers_in_parcel_data <- 
  parcel_sf %>%
  group_by(parcelnum) %>% 
  mutate(n = n()) %>%
  ungroup() %>% 
  filter(n > 1) %>%
  select(parcelnum, address, legaldesc, row_num)
duplicate_parcel_numbers_in_parcel_data

Simple feature collection with 78 features and 4 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -83.17034 ymin: 42.32219 xmax: -83.0055 ymax: 42.4416
epsg (SRID):    4326
proj4string:    +proj=longlat +ellps=WGS84 +no_defs

Plotting the parcel data shows that, within groups of parcels that have the same parcel number, some of the pacrcels are geometrically identical while others are disjoint. We thus a apply a spatial join, using sf::st_join, to find, among each of these groups of parcels, the parcels with different parcel numbers but which are geometrically the same.

#join the result with itself, on the basis of sameness of identity of spatial identity 
spatial_repeats <- st_join(duplicate_parcel_numbers_in_parcel_data,
                           duplicate_parcel_numbers_in_parcel_data, 
                           st_equals, left = FALSE) %>% filter(row_num.x != row_num.y)
#select one element from each group with the sf objects for which the polygon covers the same area
selection_vector <- 1:nrow(spatial_repeats)
for (index in 1:nrow(spatial_repeats)) {
  selection_vector[index] <- !(spatial_repeats$row_num.x[index] %in% spatial_repeats$row_num.y[1:index]) 
}
selection_vector <- as.logical(selection_vector)
#the unique parcels, as described in unique_parcels
unique_parcels <- spatial_repeats[selection_vector,]
#the unique parcels, as described in parcel_sf, with a row number added to the end of the parcel number (character vector)
unique_parcels <- parcel_sf %>% filter(row_num %in% unique_parcels$row_num.x) %>%
  mutate(parcelnum = str_c(parcelnum, "_", row_number()))
#cut out the set of sf objects that have spatial repeats  
parcel_sf <- parcel_sf %>% filter(!(row_num %in% spatial_repeats$row_num.x))
#bind the the set of unique spatial objects to parcel_sf
parcel_sf <- parcel_sf %>% rbind(unique_parcels)
#check for overlap of the parcels in parcel_sf
spatial_repeats <- st_join(parcel_sf, parcel_sf, 
        st_overlaps, left = FALSE) %>% filter(row_num.x != row_num.y)

although coordinates are longitude/latitude, st_overlaps assumes that they are planar

#create vectors or row numbers in parcel_sf whose corresponding elements specify objects in parcel_sf with overlapping geometries.
x_rows <- spatial_repeats$row_num.x
y_rows <- spatial_repeats$row_num.y
#examine a random selection of these cases
set.seed(55); sample <- sample(1:nrow(spatial_repeats), 20)
parcel_selection <- spatial_repeats[sample,] %>% select(row_num.x, row_num.y)
plots <- map(1:20, function(index){
  row_1 <- (parcel_sf %>% select(parcelnum))[parcel_selection$row_num.x[index],]
  row_2 <- (parcel_sf %>% select(parcelnum))[parcel_selection$row_num.y[index],]
  plot(rbind(row_1, row_2))
})

#remove the row_num variable from pracel_sf
parcel_sf <- parcel_sf %>% select(-row_num)
rm(duplicate_parcel_numbers_in_parcel_data, selection_vector, unique_parcels)

The dismantle permits data also contains some duplicate parcel numbers over rows, in some cases over rows that contain address information suggesting that the location is different. (It also indicates that a few individual locatations, identified with addresses, had more than one associated dismantle permit. This need not be problematic—a permit could expire before the work is carried out, or there could be more than one structure on a parcel.) These repetitions are in duplicate_parcel_numbers_over_distinct_addresses, which contains 260 rows. I should also note that, in most of these cases with duplicate parcel numbers (and addresses indicating different locations), the recorded latitude and longitude are identical. We can thus infer that some of this location information is incorrect.

#repeated parcel numbers in the dismantle permits data
dup_par_num_in_dismantle_data <-
  dismantle_permits_sf %>%
  group_by(`Parcel Number`) %>% 
  mutate(n = n()) %>%
  ungroup() %>% 
  filter(n > 1) %>%
  select(`Parcel Number`, `Site Address`) %>%
  arrange(`Parcel Number`)
rm(dup_par_num_in_dismantle_data)
#parcel numbers in the dismantle permits data that are distributed over disinct address strings
dup_par_num_over_distinct_addresses <- 
  dismantle_permits_sf %>%
  group_by(`Parcel Number`) %>%
  mutate(parcel_number_occurances = n()) %>%
  ungroup() %>%
  filter(parcel_number_occurances > 1) %>%
  group_by(`Parcel Number`, `Site Address`) %>%
  mutate(m = n()) %>%
  filter(m < parcel_number_occurances) %>%
  arrange(`Parcel Number`)
#view selected portions of dup_par_num_over_distinct_addresses
as.data.frame(dup_par_num_over_distinct_addresses) %>%
  select(`Parcel Number`, `Site Address`, parcel_number_occurances, m)

#remove extraneous variables and then geocode the dismantle site addresses
geocoded_duplicates <- dup_par_num_over_distinct_addresses %>%
  as.data.frame %>%
  select(-parcel_number_occurances, -m, -geometry) %>%
  mutate(location = map(str_c(`Site Address`, ", Detroit, Michigan"), geocode_pause))
  
#avoid having to do the geocoding again when re-running the code
write_rds(geocoded_duplicates, "./data/geocoded_duplicates.rds")  

Having identified some potentially problematic data points, we need to investigate a few of them. As in the following example.

temp1 <- read_rds("./data/geocoded_duplicates.rds")
#Isolate some variables for viewing
temp2 <- temp1 %>% select(`Parcel Number`, location, `Permit Location`)
#eleven seeminly separate locations with the same listed parcel number
temp3 <- temp1 %>% filter(`Parcel Number` == "13006809.")
#put into an sf form
temp4 <- st_as_sf(unnest(temp3), coords = c("lon", "lat"), crs = st_crs(parcel_sf))
#function for putting a set of of sf point on a satelite map
sf_points_map <- function(sf_df) {
  df <- sf_df %>% mutate(longitude = st_coordinates(sf_df)[,1],
                                latitude = st_coordinates(sf_df)[,2]) %>%
  as.data.frame %>% select(longitude, latitude)
  
  #left/bottom/right/top for bounding box
  bounding_box <- c(min(df$longitude) - 0.002, min(df$latitude) - 0.002, 
                    max(df$longitude) + 0.002, max(df$latitude) + 0.002)
  detroit_gg <- get_map(location = bounding_box, 
                        maptype = "satellite")
  ggmap(detroit_gg) + geom_point(data = df, aes(x = longitude, y = latitude)) 
}
#map of points with Parcel Number listed as "13006809."
sf_points_map(temp4)

converting bounding box to center/zoom specification. (experimental)
Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=42.410747,-83.043346&zoom=17&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false

#It is also notable that "13006809." is not reflected as a parcel number in parcel_sf (noting that parcel_sf$parcelnum is an integer vector)
parcel_sf %>% filter(parcelnum == as.integer("13006809."))

Simple feature collection with 0 features and 37 fields
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
epsg (SRID):    4326
proj4string:    +proj=longlat +ellps=WGS84 +no_defs
 [1] sev        building_s landmap    pre        taxable_st last_sale  last_sa_01
 [8] objectid   num_buildi floor_area owner1     parcelnum  sqft       owner_coun
[15] land_value legaldesc  related_pa zoning     nez        address    owner2    
[22] total_acre taxable_va owner_city depth      year_built last_terms improved_v
[29] owner_stat taxpayer   council_di owner_zip  zip_code   property_c frontage  
[36] owner_stre ward       geometry  
<0 rows> (or 0-length row.names)

#try a spacial join to see where these points in the dismantled dataset hook up with the parcels dataset
temp5 <- st_join(parcel_sf, temp4, left = FALSE, st_contains)

although coordinates are longitude/latitude, st_contains assumes that they are planar

temp5$`Parcel Number`

[1] "13006809." "13006809." "13006809."

#three rows, all with `Parcel Number` = "13006809."
#plot with crs = 3857 to match google maps
sf_points_map(temp4) +  
  geom_sf(data = temp5[1,] %>% select(geometry), crs = 3857, inherit.aes = FALSE)

converting bounding box to center/zoom specification. (experimental)
Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=42.410747,-83.043346&zoom=17&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false
Coordinate system already present. Adding new coordinate system, which will replace the existing one.

rm(temp1, temp2, temp3, temp4, temp5)

We thus see that at least some of the parcels contain multiple structures. It is also notable that some of the positions associated with the dismantle permits, particularly those in the streets, are far enough from the the parcel to potentially be closer to other parcels. In particular, investigation of the surrounding area with google earth suggests that the point in temp4 that is the furthest south is quite close to a dismntled building that is well-outside of the parcel.

We resume the development our our model.

#read in the saved set of examples that involve repeats of parcel numbers over different addresses.
dismantle_duplicates_geocoded <- read_rds("./data/geocoded_duplicates.rds")
#unpack the column of data frames (the outputs of geocode())
dismantle_duplicates_geocoded <- remove_null_locations(dismantle_duplicates_geocoded) %>%
   unnest(location)
#note the 16 rows of dismantle_duplicates_geocoded for which another row lists the same longitude and latitude
dismantle_duplicates_geocoded %>% 
  group_by(lon, lat) %>% mutate(n = n()) %>% filter(n > 1)

  
#keep only the first of any group of rows with the same longitude and latitude
dismantle_duplicates_geocoded <-
  dismantle_duplicates_geocoded %>% distinct(lon, lat, .keep_all = TRUE)

Using Google Street View to investigate some more dismantle permit entries for which there are other entries with the same parcel number but different addresses, we find mostly different locations at the same buidling (e.g. the two sides of a duplex), or perhaps two attached buildings.

Continuing our construction:

#manually cut out the one remaining apparent duplicate location
dismantle_duplicates_geocoded <- dismantle_duplicates_geocoded %>% 
  filter(!(`Site Address` == "3200 E LAFAYETTE-MARTIN LUTHER KING HIGH"))

dismantle_duplicates_geocoded <- dismantle_duplicates_geocoded %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

#It is odd that this step is necessary, given that CRSs of the two data frames appear to be the same
st_crs(dismantle_duplicates_geocoded) <- st_crs(parcel_sf)

#where possible, find the parcels that contain the positions corresponding to the coordinates
spacial_join_within <- 
  st_join(dismantle_duplicates_geocoded, 
          parcel_sf %>% select(parcelnum), 
          st_within)

#change the parcel numbers (for the dismantle permits data) to those for the parcels that contain the parcels specified by the coordinates. remove the parcelnum column (from the parcel_sf data frame)
dismantle_duplicates_geocoded <- spacial_join_within %>%
  mutate(`Parcel Number` = ifelse(!is.na(parcelnum),
                                  as.character(parcelnum),
                                  `Parcel Number`)) %>%
  select(-parcelnum)

#remove from dismantle_permits_sf all of the rows with repeated parcel numbers and then
#replace then, when possible, with geocoded rows. first make the CRSs consistent (with the observation, again, that they would appear to already be the so).
st_crs(dismantle_duplicates_geocoded) <- st_crs(dismantle_permits_sf)

#PERHAPS CHANGE THIS...
dismantle_permits_sf <- dismantle_permits_sf %>% 
  filter(!(`Parcel Number` %in% dup_par_num_over_distinct_addresses$`Parcel Number`)) %>%
  rbind(dismantle_duplicates_geocoded)

rm(dismantle_duplicates_geocoded, dup_par_num_over_distinct_addresses)

We now begin the process of creating a set of labels over the buildings. The first step will be to create a dataframe in which each row is a “building”—the set parcels (elements of parcel_sf) such that each either (1) has a building on it, as indicated in parcel_sf, or (2) a dismantle permit for the same location as the parcel. I will begin with the identification of the buildings.

buildings <- parcel_sf #%>% 

typeof(parcel_sf$parcelnum)
typeof(dismantle_permits_sf$`Parcel Number`)

parcel_join_sf_parcel_first <- 
  inner_join(parcel_sf %>% select(parcelnum) %>% mutate(parcelnum = as.character((parcelnum))),
            dismantle_permits_sf %>% as.data.frame() %>% select(`Parcel Number`),
            c("parcelnum" = "Parcel Number")) %>%
  arrange(parcelnum)

parcel_join_sf_dismantle_first <-
  full_join(dismantle_permits_sf %>% select(`Parcel Number`), 
            parcel_sf %>% as.data.frame() %>% select(parcelnum),
            c("Parcel Number" = "parcelnum")) %>%
  arrange(`Parcel Number`)

st_crs(parcel_sf)
st_crs(dismantle_permits_sf)

#Both of these yield an EPSG of 4326. Nethereless, the follow step seems necessary.

st_crs(dismantle_permits_sf) <- st_crs(parcel_sf)

temp <- st_join(dismantle_permits_sf, parcel_sf, join = st_within)

temp2 <- inner_join(temp, as.data.frame(parcel_join_sf_parcel_first), by = "parcelnum")
  
temp2 <- temp %>% filter(`Parcel Number` != parcelnum)

temp2 <- temp2 %>% select(`Parcel Number`, parcelnum, `Site Address`, address)


#note that objectid in parcel_sf has no NA values, and that objectid is a key in parcel_sf
parcel_sf %>% as.data.frame() %>% filter(is.na(objectid))
parcel_sf %>% count(objectid) %>% filter(n > 1)

#note that `Permit Number` is a also a key, for dismantle_permits_sf
dismantle_permits_sf %>% as.data.frame() %>% filter(is.na(`Permit Number`))
dismantle_permits_sf %>% count(`Permit Number`) %>% filter(n > 1)

#with the above in mind, we examine the NA values for these columns in our full join, beginning with objectid in parcel_sf. The idea is that any NA value for objectid in parcel_join would correspond to a parcel in the dismantle permits dataset that is either not reflected in the parcels dataset or for which there is a flaw in the data.
parcel_join %>% filter(is.na(objectid))

dismantle_permits_sf %>% filter(`Parcel Number` == "22108512.")

#dismantled_structures <- initial_join %>% filter(!is.na(`Permit Number`))

parcel_sf %>% filter(parcelnum == "22108512.")

parcel_tmp <- blight_violations_sf[c(1, 2),]

st_distance(parcel_tmp)

m_list <- mapply(st_distance, blight_violations_sf$geometry[1], fred$geometry)

dismantle_permits_sf %>% as.data.frame() %>% filter(is.na(`Parcel Number`))

temp4 <- st_distance(parcel_sf[195099,], 
            blight_violations_sf[22233:22255,])

st_crs(blight_violations_sf) <- st_crs(parcel_sf)

plot(blight_violations_sf, max.plot = 1)

Maps:

library(ggmap)

ggplot(blight_violations_sf %>% select(geometry)) +
  geom_sf() + 
  geom_density2d(data = blight_violations_sf %>% select(geometry))


detroit_gg <- get_map("Detroit", maptype = "toner-lite",
                      zoom = 11)

temp <- get_map()

ggmap(detroit_gg)